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Abstract 

The late-time behavior of a scalar field on fixed Kerr background is ex- 
amined in a numerical framework incorporating the techniques of conformal 
compactification and hyperbolic initial value formulation. The applied code is 
l+(l-l-2) as it is based on the use of the spectral method in the angular direc- 
tions while in the time-radial section fourth order finite differencing, along with 
the method of lines, is applied. The evolution of various types of stationary and 
non-stationary pure multipole initial states are investigated. The asymptotic 
decay rates are determined not only in the domain of outer communication but 
along the event horizon and at future null infinity as well. The decay rates 
are found to be different for stationary and non-stationary initial data, and 
they also depend on the fall off properties of the initial data toward future 
null infinity. The energy and angular momentum transfers are found to show 
significantly different behavior in the initial phase of the time evolution. The 
quasinormal ringing phase and the tail phase are also investigated. In the tail 
phase, the decay exponents for the energy and angular momentum losses at J^"*" 
are found to be smaller than at the horizon which is in accordance with the be- 
havior of the field itself and it means that at late times the energy and angular 
momentum falling into the black hole become negligible in comparison with the 
energy and angular momentum radiated toward . The energy and angular 
momentum balances are used as additional verifications of the reliability of our 
numerical method. 
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1 Introduction 



The late-time behavior of hnear fields propagating on Schwarzschild background has 
been studied in much detail since the pioneering work of Price [1] . By now it is well 
understood and there is also compelling numerical evidence supporting the theoretical 
picture. The analysis is considerably simplified by the fact that, due to spherical sym- 
metry of Schwarzschild spacetime, the perturbing field can be uniquely represented as 
a linear combination of spherical harmonics and the coefficients of the separate mul- 
tipole components evolve independently. It has been justified by several authors that 
for compactly supported, generic, i.e. non-stationary initial data, an (/,m) multipole 
component decays at late times as t"(2'+3)^ regardless of the type — scalar, electromag- 
netic, gravitational, etc — of the perturbing field. 

Since the Kerr spacetime is not spherically symmetric, the late-time decay on 
Kerr background is more complicated and the numerical results had also remained 
controversial until quite recently. Due to the lack of spherical symmetry, even if ini- 
tially only a single multipole component is excited, inevitably additional multipole 
components are also generated during the time evolution. Fortunately, the axial sym- 
metry along with the equatorial refiection symmetry of the Kerr spacetime restricts 
the allowed multipole components. Nevertheless, somewhat surprisingly, it was found 
that in certain coordinate systems, including the Boyer-Lindquist and ingoing Kerr 
ones, the decay rates of the excited multipole components show a dependence on the 
multipole indices of the initially excited mode that cannot be accounted for simply 
by referring to the symmetries of the background. For scalar perturbations, in the 
case of non-stationary initial data specifications, this sort of behavior was predicted 
by the analytic investigations of Barack and Ori [21 |3], and of Hod [1]. In addi- 
tion, predictions for the decay exponents — relevant for the case with Boyer-Lindquist 
coordinates — were also derived in [21E1I1]- Numerical justifications of these analytic 
results were also provided by Burko and Khanna OE], Tiglio, Kidder and Teukolsky 
[7], Gleiser, Price and Pullin [8], Krivan, Laguna and Papadopoulos [9], and Krivan 

m- 

To understand the main difficulties showing up in the study of the evolution of a 
linear field, $, satisfying the homogeneous wave equation 

V^V'^$ = 0, (1) 

on Kerr black hole background, note that while due to spherical symmetry of the 
Schwarzschild spacetime the family of two-spheres invariant under the action of the 
rotation group provides a geometrically distinguished and thereby unique foliation of 
the spacetime, such a unique foliation does not exist in the Kerr case. Although some 
of the most frequently used coordinates, such as the Boyer-Lindquist or Kerr-Schild 
coordinates, offer immediate choices, the pertinent topological two-spheres do not 
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coincide and, in turn, a pure multipole state with respect to one of these fohations is 
not a pure one with respect to the other, and vice versa. For a related discussion see 
e.g. [H], where the exphcit fohation dependence of the decay rates was also clearly 
demonstrated. Tiglio, Kidder and Teukolsky [7] provided numerical confirmation of 
these analytic indications by comparing evolutions based on Boyer-Lindquist or Kerr- 
Schild coordinates. Numerical investigations using Kerr-Schild coordinates were also 
done by Scheel et al. [12], with results that are in agreement with those of [7]. It was 
found, in particular, that in Kerr-Schild coordinates the late-time decay rates relevant 
for the lowest admitted multipole components do depend on the multipole indices of 
the exciting pure initial data only up to the extent that can be deduced from the 
axial and reflection symmetries of the Kerr spacetime. Note, however, that even in 
this simple case the decay rates of the higher multipole components differ from the 
corresponding decay rates in the Schwarzschild case. In addition to the non-unique 
choice of the two-spheres, the consequences of using different time slicings and the 
associated loss of the generality of the initial data specifications also increase the 
complexity of the picture. A significant contribution to the clear up is due to Tiglio, 
Kidder and Teukolsky [7] who carried out a systematic comparison of the decay rates 
based on the use of Boyer-Lindquist and Kerr-Schild coordinates. 

In favoring the non-Kerr-Schild type coordinate systems, Burko and Khanna [6] 
provided numerical justification of the fact that the decay rates of tails in Boyer- 
Lindquist and in ingoing Kerr coordinates agree, in spite of the fact that the time 
slices are different as the latter are horizon penetrating ones. They also claimed 
that these coordinate systems belong, in a certain sense, to an equivalence class in 
which pure multipole initial data are distinguished and their indices are uniquely 
determined by the late time decay rates yielded by them. By making use of the 
relation between the Boyer-Lindquist and ingoing Kerr coordinates — see equations 
© - © — they argued in [6] that a spherical harmonic mode with indices (/, m) in 
Boyer-Lindquist coordinates is also a spherical harmonic mode with the same values 
for {I, m) in ingoing Kerr coordinates, which may explain the similar behavior of 
the tails in these two coordinate systems. The latter point, however, needs further 
justification and more rigorous arguments. 

Providing a verification of the proposal of Burko and Khanna O [6] , in a slightly 
different setup, Zenginoglu and Tiglio |13] carried out a numerical investigation of the 
decay rates based on the use of conformal compactification, such that the compactified 
versions of both the Boyer-Lindquist and the ingoing Kerr coordinates were applied. 
Although their results for the asymptotic decay rates are not very extensive, the 
reported values are in good agreement with that of [H [5l [6]. In particular, they 
support the suggestions of Burko and Khanna [6]. One of the main advantages of 
using conformal compactification is that the decay rates at future null infinity, J^^, 
may also be determined. The pertinent decay rates reported in [13] are in agreement 
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with the predictions of [1]. 

The aim of this paper is to strengthen the developing consensus by presenting 
the results of our numerical investigations on the late-time evolution of tails in the 
Kerr spacetime. The applied analytic setup is based on the use of the techniques of 
conformal compactification and hyperboloidal initial value formulation. To simplify 
the presentation and also to provide a straightforward way to derive the basic equation 
to be solved — instead of a technical but covariant argument — we have chosen a non- 
covariant one in which two subsequent coordinate transformations are applied. In 
addition to the usually chosen quickly decaying initial data specifications some slowly 
decaying cases, which have not been studied yet, with fall off for $ toward 
are also considered. Similarly, to get some more insight about the time evolution of 
the scalar field on Kerr background, not only the axially symmetric case but in a few 
cases excitations with azimuthal indices m = ±1 are also studied. In addition to the 
usual convergence criteria the energy and angular momentum conservation laws are 
also monitored. Our results concerning the associated energy and angular momentum 
transports also provide some important clue about the evolution of linear fields on 
rotating black hole backgrounds. 

The paper is organized as follows. In Section the applied analytic and numerical 
framework is introduced. In Section [3l our main results are presented, and Section |4] 
contains our final remarks. Throughout this paper, geometric units, with c = G = 1, 
are applied. 

2 Preliminaries 

This section is to recall some of the basics of the applied mathematical and numerical 
framework. 

2.1 The analytic setup 

The Kerr metric in Boyer-Lindquist coordinates (t, r, 9, 0) reads as 

dt^ ^ dtdcf) + —dr^ + Srf^^ + - sin^ ^d</)^ (2) 

where E = + ^os^ 6, A = r"^ - 2Mr + a^, r = (r^ + a^f - a^A sin^ 6. M and 
a denote the mass and the angular momentum per unit mass of the Kerr black hole, 
respectively. 

The desired new coordinates, which have the advantages that the time slices are 
horizon penetrating and connect to future null infinity such that no boundary condi- 
tions are needed, are introduced by performing the following two subsequent coordi- 
nate transformations. 



, 2 / 2Mr\ 
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(1) First, the time and azimuthal coordinates, t and 0, of the Boyer-Lindquist 
system are replaced by r and ip defined as 



r = t-r + J dr — - — (3) 
^ = 4> + J dr^. (4) 

The coordinates (r, r, 6, (f) are usually referred to as ingoing Kerr coordinates. 
As opposed to the Boyer-Lindquist coordinates, these coordinates remain reg- 
ular at the outer horizon "H^ separating the domain of outer communication 
from the interior domain between the inner and outer event horizon. These re- 
gions are indicated by domains I and II in Figure [H respectively. As the radial 
coordinate, r, is not changed, the inner and outer event horizons, H- and 
are located at r = = M =p a/ — a^, respectively. As indicated in Figure 
[H the r = const time-slices intersect 71+, justifying the terms 'ingoing Kerr 
coordinates' and 'horizon-penetrating slicing'. Note also that in the r — oo 
limit the r = const time-slices tend to spacelike infinity, denoted by 

(2) Second, we replace the ingoing Kerr coordinates r and r by the new time coor- 
dinate T and by the compactified radial coordinate R, defined via the implicit 
relations 

T = T+ii|^-4Mlog(|l-iJ'|) (5) 

(6) 



The transformation ([5]) - ([6]) yielding the new coordinates (T, R, 6, ip) differs 
from the conformal transformation proposed by Moncrief [H] (see also [T^ 
[TB] ) only in the logarithmic term in which — according to the conclusions 
of [T7] — is inevitably present in most of the analogous conformal scri fixing 
compactifications of the Kerr spacetime. The T = const surfaces (see Figure 
[1] for an illustration), tend to future null infinity (J^^) as r — )■ oo, which is 
represented by the R = 1 boundary in the new coordinate^ and are spacelike 
in a neighborhood of the domain of outer communications. It is straightforward 
to see that the coordinate basis fields {d/dT)^ and [d/dip)^ are Killing vectors. 
A major advantage of these coordinates is that one can place the inner and outer 



^Note that some of the coordinate components of the metric tensor blow up at i? = 1, which 
in the framework of conformal compactification is usually compensated by replacing the physical 
metric g^i, by a non-physical one g^y = fl^g^^, where fl is the conformal factor. In the particular 
case considered here, Q — 
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boundaries of the domain of numerical investigations inside the black hole region 
and to future null infinity. Since l-L^ and J^^ are null hypersurfaces admitting 
the waves to leave the domain of outer communication but not allowing any 
to enter, such a choice considerably reduces spurious boundary effects in the 
domain of outer communication (see also [18] for a discussion on the effects 
of boundary conditions). Furthermore, there is no need to impose any sort of 
sophisticated boundary condition and we find it to be completely satisfactory 
to require simply the field equation to hold at the edges. In our numerical 
calculations, the inner boundary (see the thick dotted line in Figure [1]) is placed 
slightly beyond the event horizon and the outer boundary precisely at i? = 1. 

In addition to the above transformations, we also introduce the rescaled field 
variable 

The pertinent form of ([T]) for \1/ is regular in both of the spacetime regions I and II, 
it remains reg ulail up to and including the boundary at i? = 1, and it can be given 
in the form 

1 



+ ar i9r + an Or + a^d^ + aQ + qa^s^) ^ 



where Y2 = ¥2(0, (f) = iy'^(3cos^0 — 1) is the spherical harmonic function with 
/ = 2, m = and A52 stands for the spherical Laplace operator 

1 d f d\ 1 d"^ 
sm 6 06 \ 06 J sm 6 Otp^ 

It is an important property of ([HD that all the coefficients ajm, qtr, ot^, o/j^, qt, 
an, a^, oq, oa, ajr ^tt functions of R exclusively, i.e. they do not depend 
on T, 6 and (f. Explicit expressions for these coefficients are given in Appendix A. 
From equation (D.14) of |2S] and from the vanishing of the scalar curvature of the 
Kerr metric, it follows that gf^'V/y^^ = VL^{g^"'W^Vu - \S)'^ , where V^, is the 
covariant derivative corresponding to the rescaled metric g^y = Vt^g^y, S is the scalar 

^It is known that a conformal compactification, which is achieved in the present case by the 
introduction of the coordinates (T, R, 6, ip), also needs to be compensated by a rescahng of the type 
([7]) in the field variable, otherwise the relevant form of ([ij for $, exactly as the physical metric 
does, becomes singular at (See also [TH] for an interesting related argument applicable in 

investigations carried out in the frequency domain.) 
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Figure 1: At the left lower block a schematic version of the Carter- Penrose diagram 
of maximal analytic extension of the Kerr spacetime is shown while on the top a part 
of this diagram is enlarged to show in more detail regions I and II. In these regions, 
the Boyer-Lindquist t and r coordinate lines are indicated by the thin dotted lines. 
Some of the r = const and T = const time-slices are indicated, for the case with 
M = 1 and a = 0.5, by the solid blue and dashed red lines, respectively. The red 
lines illustrate that the (T, R) coordinates are horizon penetrating and connect to 
J^"*". It is worth mentioning that constant time-slices of the Kerr-Schild coordinates 
coincide with that of the ingoing Kerr coordinates and these two coordinate systems 
differ only in the applied foliation by topological two-spheres. 
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curvature of g^^^, g^^ = VL~'^g^^ and Vt^^ = -jz^- Using this equation, the relation 

ci'O = 6^TT'5'(a^^ + aip^Y2) can be obtained. Note that while depends only on 

R, S also depends on 6. The reason for writing the field equation in the form (|8]), 
where it is solved for (9tt^, is that it is this form that is needed for the application 
of the Runge-Kutta method, which is the method that is used to compute the time 
evolution of the field. It is a further important property of ([H]) that the denominator 
on the right hand side is non-vanishing in a neighborhood of the domain of outer 
communication. 

2.2 The numerical method 

In solving ([H]), the angular dependence of the field is handled by spectral decomposi- 
tion based on the foliation by two-spheres given by the t,r = const — or equivalently, 
T,R = const — surfaces. Accordingly, \1/ is replaced by the series of spherical harmon- 
ics 

^ = j2i'riR,T)Yr{9,^), (10) 

which reduces the evolution problem to determining that of the coefficients ipf^{R, T) 
in the T — R section. In practice, the series ffTOj) has to be truncated to a finite sum; 
nevertheless, as it converges quickly, the desired precision can always be guaranteed 
by keeping sufficiently many terms. Note that (jSj) has a relatively simple structure, in 
which only the denominator, in particular its second term, is responsible for a mixing 
between different spherical harmonic modes. Following the basic ideas described in 
more detail in [201 [21] , the division by the term -|- a^/^ Y2 is converted into a 
multiplication by applying the identity 

1/(1 + x) = l-x + x^ -x^ + -x^ ■■■ , (11) 

along with the choice 

^2 I'^TT °f relation 

On the right hand side of ffTTj) . one has to keep as many terms as needed to get the 
error induced by this approximation to be smaller than a fixed tolerance. In our 
simulations the error tolerance is chosen to be 10"^", corresponding to the extended 
precision arithmetic that we use. For a = 0.5, which is the value that is taken in most 
of our calculations, 11 terms are kept in (fTTj) . which is sufficient since la^-^Fg^/'^TTl < 
0.00775 holds for a = 0.5. The coefficients G^/^Tg'™', appearing in are closely 
related to Clebsch-Gordan coefficients and are described in standard mathematical 
references. In addition to the above replacements the well-known relations 

= -/(/ + l)y;™ and d^Y/^ = \mY{^ (13) 
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are also applied. Due to the axial symmetry of the Kerr background and the linearity 
of dH]), the time evolution of the coefficients ipplR, T) with different azimuthal indices 
decouple, which, in particular, allows one to investigate their evolution separately. 
As the background also possesses the equatorial reflection symmetry, the coefficients 
ipp{R,T) with odd differences in their / indices also decouple, which considerably 
simplifles the study of the evolution of pure multipole initial data. In our simulations 
pure multipole initial data with /' < 5 are considered, where /' is used to denote the 
polar index of the initial data. It is found that with such initial data (the complete 
speciflcation of our initial data can be found in Section [3?T]) it is sufficient for the study 
of several modes up to / = 4 to keep only the first 8 non-zero terms in the multipole 
expansion ffTOl) of due to the fast convergence of this series. This means that the 
multipole components taken into account are those for which I = \m\, . . . , \m\ + ((/' — 
|m|) mod 2) + 14, since / > |m| for any spherical harmonic function and the non-zero 
modes are those for which /' — / is even. 

(2) 

As the coefficient Qjrprp IS proportional to the angular momentum per unit mass 
parameter of the Kerr spacetime, it follows immediately from the above discussion 
that the spherical modes completely decouple in the Schwarzschild limit. It is also 
worth keeping in mind that, in virtue of ([3]) - dS]), a spherical harmonic mode with 
indices (/, m) defined with respect of the Boyer-Lindquist coordinates is also a spher- 
ical harmonic mode, with the same (/, m) indices, with respect to the compactified 
ingoing Kerr coordinates (T, i?, 9, ip), and vice versa. 

In solving ([8]), it is recast into a first order form by introducing the new auxiliary 
variables = Ot"^ and = (9/j\l/, and in the T — R section the method of lines 
is applied to compute the time evolution of the multipole coefficients of and 
Fourth order finite differencing is used in the radial direction while for the 
time evolution a fourth order Runge-Kutta scheme is applied. In order to suppress 
high frequency instabilities, a standard fifth order Kreiss-Oliger dissipation term, as 
proposed in [22], is added to the right hand side of the equations, and the constraint 
= dfi'^ is imposed after each time step. It is found that 2048 grid cells are 
sufficient in most of our simulations, although in a few cases, e.g. in checking the 
convergence, grids of size 1024, 4096 and 8192 are also used. The time step is chosen 
to be ten times the radial lattice spacing in accordance with the smallness of the 
coordinate light speed. 

2.3 Conserved currents 

As {d/dT)'^ and [d/dipY are Killing vectors, the energy and angular momentum 
currents 




and 




(14) 
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are divergence free, where 

%u = (V^<I>)(V,$) + (V.$)(V,,$) - (7^,(Va$)(V^$) (15) 

is the energy-momentum tensor. It follows then from Stokes' theorem that for a 
spacetime domain N with boundary c^A^ and outward pointing unit normal vector 
at dN , the energy and angular momentum balance relations 

f n^E''= I ^^£" = and [ n^M^ = [ V^M'^ = (16) 
JdN Jn Jon Jn 

hold. We use these relations — in addition to the usual convergence checks — to verify 
the correctness and reliability of our code. In evaluating the integrals, is always 
chosen to be a cylindrical domain, determined by the type of relations Ti < T < T2, 
Ri < R < R2. Further details on the calculation of the integrals Jgi^n^Ef^ and 
J^j^n^M'^ are given in Appendix B. 



3 Results 



3.1 Initial conditions 

In our numerical investigations the initial data is always 'pure' in the sense that only 
one of the coefficients ipj^ or its time derivative is chosen to be non-zero. We consider 
eight different types of initial data labeled by numbers 1, . . . , 8. The type 1 and 2 
initial data are compactly supported whereas for the others the physical field $ have 
various fall off property toward future null infinity. 

Type 1 initial data is such that \i/ = and, for specific values of the indices (/, m), 
dT'4'r is chosen to be times the 'bump' function 

^ ^ jexp (-p^z^ - + ^) , ifc-w/2<R<c + w/2, ^^^^ 

[ , otherwise , 

where w denotes the width of the bump and c determines its center. ^ is smooth 
and it is of compact support. We refer to this initial data as non-stationary as the 
time derivative of the field is non-zero. 

Type 2 initial data are complementary to type 1 in the sense that ipp is chosen 
to be times the bump function determined above, while Ot"^ is chosen to be 
identically zero. This type of initial data is going to be referred as stationary since 
the time derivative of the field vanishes. 

For type 1 and 2 initial data, in all of our numerical calculations, the values w and 
c are chosen to be 0.1 and 0.7, respectively, i.e. the bump is narrow and it is located 
close but outside the outer event horizon. 
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The non-compactly supported initial data specifications, type 3-8, also comprise 
complementary pairs. Namely, type 3, 5 and 7 are non-stationary while type 4, 6 and 
8 are the corresponding stationary ones yielded from 3, 5 and 7 by interchanging the 
initial data specifications following the rule ^1/3,5,7 — 57^^4,6,8 and (9r^'3,5,7 — )■ ^^4,6,8- 

Type 3 and 4 initial data are specified by setting either the coefficient i/jJ^ or its 
time derivative to be identically one. In virtue of ([7]), ^ = r $, i.e. in this case $ or 
its time derivative possesses 1/r fall off toward J^~^. 

In the construction of type 5-8 initial data, the bump function ^ with parameters 
c = 0.9 and w = 0.4 is applied, which vanishes for R < 0.7, i.e. near to the outer event 
horizon, and it is non-zero at i? = 1. For type 5 and 6 either the coefficient ip]^ or its 
time derivative is chosen to be 1 — i? times ^] thereby, it vanishes linearly at i? = 1 
while it is zero in a neighborhood of the event horizon. Type 7 and 8 initial data 
specifications are constructed analogously multiplying with parameters c = 0.9 
and w = 0.4, by the factor (1 — R)"^. 

It is important to note that type 5, 6 and 7, 8 initial data specifications imply 
and fall off property toward null infinity for $ or 9^$, respectively. In the 
Schwarzschild case, initial data with the power law fall off property of type 5, 6 or 7, 8 
initial data were studied by Price [T] — referred by him as 'static' initial data — and 
Leaver |23]. Note also that in [H [23] the pertinent decay rates were found to be 
smaller than the corresponding rates for initial data with compact support. 



3.2 Power law decay exponents 

After an initial burst phase, the field goes through a longer lasting quasi-normal 
ringing period while the evolution is closed by a power law decay. This, in particular, 
means that for large values of T the multipole expansion coefficients, at constant R, 
decrease as ~ T~", where the value oi n = n{l', I) is found to be integer in analytic 
investigations [H El El IH [HI [2S1 [2B]- Here, using the terminology introduced in [3, [Hj, 
the primary multipole index of the pure initial data is referred to as I' while that 
of the examined excited multipole modes as /. Since the azimuthal mode m is fixed 
during the evolution, its value will not be explicitly indicated. 

The power law decay, in general, starts at T values of a few hundred, where T is 
measured in mass units. In numerical investigations the decay exponents n = n{l', I) 
are extracted by evaluating the local power index (LPI) 



dlnT 



In all of our numerical investigations the mass of the black hole M is chosen 
to be one, while in most of the cases the choice a = 0.5 is applied, although some 
simulations with a = 0.2 and a = 0.995 are also done to justify that the decay rates 
do not depend on the value of a. Since in former numerical investigations almost 
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exclusively the m = case has been considered, for the sake of comparison and for 
simplicity m = is chosen also in most of our simulations. Nevertheless, for type 1 
and 2 initial data specifications the cases m = ±1 are also investigated. 

Table [1] below shows the decay exponents that we find with type 1 (i.e. compactly 
supported, non-stationary) initial data. For each value of /, two numbers appear 
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Table 1: Decay exponents n{l', I) for type 1, non- stationary initial data with compact 
support, with / = 0, 1, 2, 3, 4 and /' = 0, 1, 2, 3, 4, 5. Here and elsewhere, in each slot 
for a specific value of /, the first number refers to the decay exponent found for R < 1 
while the second entry to that at i? = 1. The modes excluded by equatorial reflection 
symmetry are indicated by the sign while those which remain ambiguous by the 
sign " X " . 

in each entry. The first one indicates the decay exponent at constant R with R < 
1 (including thereby the outer horizon, 71+) while the second number signifies the 
decay rate at i? = 1 (representing future null infinity, J^~^). For R < 1, the decay 
exponents are measured at several values of R ranging from the inner boundary of 
the computational domain to J^~^. The inner boundary of the computational domain 
is set to R = 0.58, which is behind the event horizon that is at i? ~ 0.59864 for 
a = 0.5. The slots for those exponents that we cannot determine unambiguously — for 
more details see the related remarks below — are marked by the sign " x". The non- 
stationary and compactly supported type 1 initial data may be interpreted as data 
falling off 'extremely fast' toward null infinity. These type of data were used in most 
of the former numerical simulations and our results for R < 1 are in agreement with 
those obtained in El [71 El El [13] • In particular, the numbers in the first entries 
support the rules n = -|- / -|- 3 for / > /' and n = -|- / + 1 for / < /' proposed in [5] . 

The numbers in the second entries, indicating the decay rates at future null infinity, 
are in agreement with those of [13]. Nevertheless, it is worth mentioning that our 
results provide extension to that of [13] as there only the cases of type 1 initial data 
with /' = . . . 4 were considered and the decay rates of the entire composite field (not 
the individual multipole components) were determined. Our results do also fit to the 
rule n = Z -t- 2 for / > Z' and n = /' for Z < /' — 2, which was obtained analytically 
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in It is interesting, although it is in accordance with this latter rule, that for the 
excitations /' = 4 and /' = 5 two of the lowest allowed excited modes, with I = 0,2 
and / = 1, 3, have coinciding decay rates, n = 4 and n = 5, at J^~^, respectively. 

For initial data of type 1, the cases m = ±1 are also investigated. For these cases 
the / = mode is excluded but otherwise — in agreement with [5] — we find the same 
decay exponents to apply as listed in Table [H 

As mentioned above, there are certain cases where we cannot determine the decay 
exponents unambiguously. In the particular cases with /' = 2, / = 2, 4 and with /' = 3, 
I = 3, the reason is that the value of n is found to depend on the location, i.e. on 
the particular value of R for R < 1. In particular, for I' = 2, I = 2, with m = 0, 
the decay exponent is found to be 5 near the black hole while it takes the value 7 for 
larger R (exceeding the value 4 characterizing the decay rate at future null infinity). 
Note that the value 7 fits to the rule n = /' + / + 3, whereas 5 corresponds to the 
rule n = /' + / + 1. In contrast to the m = case, the decay exponent is found to 
be 7 regardless of the value i? < 1 for the cases m = ±1. Analogously, for = 2, 
/ = 4, we find the exponent 7 close to the black hole and 9 far away from it, and 
for m = ±1 we find that the exponent is 9 independently of the value of i? < 1. 
For /' = 3,/ = 3,m = 0, we find again the exponent 7 close to the black hole and 
9 far away from it. Interestingly, and in contrast to the Z' = 2, / = 2,4 cases, we 
find the same type of splitting of the decay rates for m = ±1. This behavior of 
the exponents is confirmed by the fact that if the grid size is raised from 2048 to 
4096, then the data improve similarly as shown in Figure Ei These findings definitely 
deserve further investigations, which, however, are beyond the scope of the present 
paper as they require the use of considerably longer time evolution and, in turn, much 
higher numerical accuracy. 

Table [2] shows the decay rates found for type 2 (stationary and compactly sup- 
ported) initial data. 
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Table 2: Decay exponents n{l', I) for type 2 compactly supported stationary initial 
data. 

For / > I' they are the same as with the first type of initial data, but for / < /' 
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they are greater by one, i.e. the rule n = I' + 1 + 2 appears to hold. Analogous results 
for Schwarzschild black holes were obtained in [221 ESI 123 EH]- The value n = 6 
for /' = 4, / = was also found in [6] in Boyer-Lindquist coordinates. It is worth 
mentioning that for type 2 initial data the cases of m = ±1 are also examined and 
the same exponents are found as for the m = case. For /' = 2, / = 2, 4, m = 0, a 
similar dependence of the decay exponents on i? < 1 is found as with the first type 
of initial data. For /' = 2, / = 2, the exponents 6 and 7 are found near the black hole 
and far from it, respectively. For Z' = 2, / = 4, the exponents 8 and 9 are found. For 
m = ±1, the exponents, also similarly as with the first type of initial data, are found 
to be independent of the value of i? < 1. Their values are 7 and 9 for / = 2 and / = 4, 
respectively. 

Figure [2] shows the local power index seen at various distances from the black hole 
for I' = 2 and / = 0. Figure [3] illustrates the convergence of the local power index 




Figure 2: The time dependence of the local power index is shown, for type 
2 initial data with I' = 2, I = 0, m = and a = 0.5, at R = 
0.58, 0.8, 0.9, 0.95, 0.97, 0.98, 0.99, 0.995, 0.997, 0.998, 1.0. The absolute value of the 
LPI decreases as R increases. The inner boundary of the computational domain is 
aX R = 0.58, while the event horizon is located at i? ~ 0.59864. The LPI at i? = 1 
approaches —3 as T increases, whereas for R < 1 the LPI approaches —4. The grid 
size applied is 8192. 



as the number of grid cells is increased. At a certain time, the discretization errors 
become large in comparison with the proper amplitude of the multipole component, 
which is a decreasing function of T, causing the deviation of the LPI from its correct 
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value. The error remains small for a longer time if the number of grid cells is increased. 
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Figure 3: The time dependence of the local power index relevant for grid sizes 1024, 
2048 and 4096 is shown for type 2 initial data with /' = 4, / = 0, m = and 
a = 0.5 at -R = 0.9, to illustrate the convergence of the local power index as the 
number of grid cells together with the time resolution is increased. Errors caused 
by the discretization, resulting here in the dropping down of the LPI for larger T, 
start to dominate later for larger numbers of grid cells. The number of non-zero 
multipole components taken into account is kept fixed at 8, corresponding to the 



cutoff 



14. 



Table |3] shows the decay rates found for type 3 non- stationary initial data. By 
comparing these decay exponents to the corresponding ones in Table [1], it can be seen 
that each of the decay exponents found for type 3 initial data is smaller by one than 
the corresponding exponent for type 1 initial data. In the Schwarzschild case, similar 
results were found in [H for slowly decaying (they are referred as static in [T]) 
initial data specifications. It is also important to note that Dafermos and Rodnianski 
[21] in their recent analytic investigations derived 1.5 as a lower bound on the decay 
rate of an axial symmetric field in the domain of outer communication, i? < 1 in our 
case, and 0.5 at i.e. at i? = 1, which are just below by one-half with respect to 
the lowest exponents 2 and 1 indicated in Table [3] for the /' = and /' = 2 cases, 
respectively. 

Table H] shows the decay rates for type 4 initial data. Notice that the entries are 
also smaller by one than those for type 2 initial data, listed in Table [21 with the only 
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Table 3: Decay exponents n(/', I) for type 3 non- stationary initial data with 1/r decay 
for $ at J^+. 
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Table 4: Decay exponents n(/', /) for type 4 stationary initial data with 1/r decay for 
$ at y+. 
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exception of the /' = case when the decay exponents turn out to be identical for 
type 2 and 4 initial data specifications. 

Tables E] - [H] show the decay rates that we find for type 5-8 initial data. 
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Table 5: Decay exponents n{l',l) for type 5 non-stationary initial data with 
decay for $ at J^~^. 
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Table 6: Decay exponents n{l',l) for type 6 stationary initial data with decay 
for $ at J^+. 

The following simple observations can be made about these results. 

(1) The decay exponents in Table [S] and [7| differ by 1 or from the exponents listed 
in Table [H and the same rule applies to Tables El El and [2J 

(2) The decay exponents listed in Table are smaller than or equal to those listed 
in Table [71 and the latter are smaller than or equal to those listed in Table [TJ 
The same rule applies to Tables [6l [8] and [21 

(3) There are special cases, e.g. /' = 4 for type 5 initial data, when the dominant, 
i.e. most slowly decaying, mode at infinity is definitely not the lowest allowed 
one. 
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Table 7: Decay exponents n{l',l) for type 7 non-stationary initial data with 
decay for $ at J^~^. 
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Table 8: Decay exponents n{l',l) for type 8 stationary initial data with decay 
for $ at J^+. 
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The second observation may be supported by the following argument. If we add 
type 1 and type 5 initial data, then the result is also type 5, in the sense that it also 
vanishes as 1 — i? at i? = 1, i.e. $ falls of as at J'^] thus, it can be expected to 
yield the decay exponents listed in Table [51 On the other hand, due to the linearity of 
the field equation, the late time decay rate yielded by the superposed initial data for a 
given multipole component is the smaller one of the decay rates yielded by the initial 
data sets that are superimposed. Thus, the decay exponents yielded by type 5 initial 
data have to be smaller than or equal to the exponents yielded by type 1 initial data. 
The same reasoning applies to the pairs of type 1, 7 and type 5, 7 initial data, and 
to the stationary initial data types as well. The assumption on which this argument 
relies is that the only difference between type 1, 5 and 7 (or 2, 6 and 8) initial data 
that is relevant for the decay rates that they yield is their rate of vanishing at i? = 1. 
This argument also implies that initial data specifications for \& that tend to zero 
as (1 — i?)*^, k > 1 — or those which for $ have the fall off property l/r'^^^ toward 
J^"*" — yield decay exponents that monotonically increase with k, and are bounded 
from above by the pertinent decay exponents listed in Table [1] or [2l 

3.3 Energy and angular momentum balance 

As mentioned above we use the energy balance relation to verify the correctness and 
reliability of our numerical code. Figure S] shows the time dependence of the quantity 



for various grid sizes. Region is chosen to be the cylindrical domain with Ti = 0, 
T2 = T, Ri = and R2 = 1, where denotes the location of the outer event 
horizon, while Eq denotes the initial energy 



associated with the part of the T = hyperboloidal time slice contained by the 
domain of outer communication. In Figure IH the values of 6E corresponding to 2048 
and 4096 grid cells are multiplied by 16 and 256, respectively. The applied initial data 
are of type 2 with /' = 4 and m = 0. The exact value of 6E is 0; thus, its numerically 
calculated value should be small and tend to zero as the grid size increases. Note that 
Figure H] does also justify that the convergence rate of our numerical simulations is 
at least four corresponding to the order of the applied finite differencing method in 
the T — R section. 

Figure [5] shows the time dependence of the total energy of the field in the domain 
of outer communication, along with the energies flowing through the event horizon 




(19) 




(20) 
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Energy conservation test 
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Figure 4: The time dependence of the value of characterizing the numerical 
violation of energy conservation, is shown for type 2 initial data with V = A and 
m = and for various grid resolutions. 5E appears to converge to at somewhat 
higher than fourth order rate. The number of non-zero multipole components taken 
into account is kept fixed at 8 as usual, corresponding to the cutoff Imax = 14. 



and through c^"*", i.e. the integral n^E^^ of the energy current density at i?+ and 
i? = 1, respectively — during the initial phase of the time evolution — for type 2 initial 
data with /' = 4 and m = 0. As it is clearly indicated by the graphs the energy of the 
initial bump splits into two roughly equal sized energy bumps. The inner one reaches 
the event horizon earlier and a little bit later the outer one leaves the domain of outer 
communication through J^"*". These two bumps carry away almost the entire of the 
excitation energy. 

FigureElis to show the time dependence of the corresponding quantities during the 
quasi-normal ringing phase. In this phase considerably larger fraction of the energy 
is radiated out to than swallowed by the black hole. The total energy of the field 
in the domain of outer communication decreases as ~ exp(— 2a;T), where uj ~ 0.094 
is the decay constant for the field. 

During the tail phase the energy flow through the event horizon and through J^^ , 
respectively, follows a power law ~ for some integer n. At the event horizon, 
with m = and /' even the relation n = 2n -|- 2 holds, where n denotes the decay 
exponent of the / = mode of the field, while for /' odd n = 2n + 1, where n denotes 
the decay exponent of the / = 1 mode. At J^"*", regardless of the parity of the 
rule n = 2n + 2 applies, where n denotes the decay exponent (at of the most 
slowly decaying mode. For m = ±1, the rule n = 2n -|- 1 is found to apply at the 
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Figure 5: The time dependence of the normahzed total energy in the domain of outer 
communication (sohd blue line), along with the energies flowing through the event 
horizon (dashed green line) and through J"^ (dash-dotted red line) integrated over 
time, during the initial part of the time evolution, for type 2 initial data with /' = 4 
and m = 0. Most of the energy travels in two relatively well-localized and roughly 
equal-sized bumps, moving toward the event horizon and toward J'^ ^ respectively. 
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Energy balance in the quasi-normal ringing phase 
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Figure 6: The time dependence of the normahzed total energy in the domain of outer 
communication (sohd blue line), along with the energies flowing through the event 
horizon (dashed green line) and through J^"*" (dash-dotted red line) integrated over 
time from 70 to T, is shown for a part of the quasi-normal ringing phase of the 
evolution for type 2 initial data with /' = 4 and m = 0. The energy in the domain of 
outer communication decreases exponentially. More energy goes out to infinity than 
falling into the black hole. 
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event horizon regardless of the parity of where n denotes the decay exponent of 
the lowest mode, and fi = 2n + 2 at J^+. We should also mention that there are a 
few cases in which we cannot determine the value of n properly. Since the relevant 
decay exponents at are smaller than those at the event horizon, these rules imply 
that in the tail phase the energy flowing into the black hole becomes, as T increases, 
negligible in comparison with the energy flowing out to . 

Before turning to the description of the angular momentum transports, note that 
as each term of the contraction ra^M^, appearing in the second relation of f|T6|l . 
is proportional to either 9<^$ or its complex conjugate, the study of the angular 
momentum conservation and transports is interesting only if m 7^ 0. 

Regarding the angular momentum conservation, very similar results to those 
shown in Figure S] are obtained. The overall features of the angular momentum 
transports in the quasi-normal ringing phase and in the power law decay phase are 
also similar to that discussed for the energy transports. In particular, during the 
latter phase, the power law decay ~ T"" of the angular momentum flow with n = 2n 
is found to apply at the event horizon and n = 2n + 2 at J^~^, where n denotes the 
decay exponent of the dominant mode of the field at the event horizon and at J^+, 
respectively. 

During the initial phase, the time dependence of the angular momentum transfers 
(see Figure [7]) is markedly different from the time dependence of the corresponding 
energy transports, which are similar to those depicted by Figure |5] in the m = case. 
To make the differences to be more transparent, in Figure [H] the angular momentum 
density as a function of R for certain time slices is shown. These graphs indicate 
that almost all of the initial angular momentum moves toward the black hole, while 
only a small amount moves outward. The outward moving angular momentum is 
made up of a positive and a negative part. As time passes, the density profile of the 
outgoing angular momentum spreads out in space, giving rise to the particular time 
dependence of the outflowing angular momentum at shown in Figure [71 In spite 
of the fact that almost all (97.3%) of the energy of the system is lost by T ~ 8, in 
virtue of Figure [7] about 30% of the initial angular momentum is still present in the 
domain of outer communication at that time. 

It should be noted that in the considered case, as indicated by Figure [TJ the to- 
tal angular momentum radiated out to infinity tends approximately to while T 
increases. However, we find that the asymptotic value, to which the angular momen- 
tum radiated out toward tends, depends for instance on the center of the initial 
bump. 
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Angular momentum balance 
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Figure 7: The time dependence of the normahzed total angular momentum in the 
domain of outer communication (solid blue line), along with the angular momentum 
flowing through the event horizon (dashed green line) and through J"^ (dash-dotted 
red line) integrated over time, during the initial part of the time evolution, for type 
2 initial data with /' = 4 and m = 1. In contrast with the energy, the angular 
momentum distribution becomes less localized as time passes. Approximately all of 
the initial angular momentum falls into the black hole in this case. 
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Figure 8: Angular momentum density as a function of -R at T = 1.025, T = 2.05, 
T = 4.1 and T = 5.84, for m = 1, /' = 4 and type 2 initial data. The outward moving 
angular momentum is made up of a positive and a negative part, the sum of which 
is approximately zero, suggesting that the outward moving part of the field is made 
up of counterrotating shells. The density profile of the outgoing angular momentum 
spreads out in space. 
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4 Conclusion 



Our numerical results, concerning the evolution of scalar fields on rotating black 
hole backgrounds, have been presented. The analytic and numerical scheme we have 
chosen is based on the techniques of conformal compactification and hyperboloidal 
initial value formulation which provided the possibility that the fall off properties of 
the scalar field along the event horizon as well as at future null infinity could also be 
investigated. 

It should be mentioned here that these techniques had already been applied in 
studying the evolution of linear fields on Kerr background in [13]. However, our 
results are more extensive not only because we used a high variety of the initial data 
specifications, but also because the applied 1 + (1 + 2) decomposition made the study 
of evolution of the individual multipole components possible, whereas in [13] only 
the evolution of the entire composite field could be investigated. Another technical 
difference is that while in [13] spacelike matching of regions with different properties 
of the applied conformal factors was done — this was admittedly a significant source of 
numerical error there — we applied a single analytic horizon penetrating hyperboloidal 
foliation which did not require the use of this type of matching. 

Our main results can be summarized as follows: 

In most of the studies of the decay rates of a scalar field on a black hole background 
investigations are almost exclusively restricted to the case of initial data which have 
extremely fast decay rates toward infinity. To extend the scope of numerical studies 
we applied, besides compactly supported initial data specifications, which clearly 
belong to the conventional ones, other type of initial data specifications possessing 
(with /c = 1, 2, 3) fall off toward future null infinity. In addition, we also studied 
for each subcases both stationary and non- stationary initial data. We would like to 
emphasize, however, that even in the conventional case, all of our results concerning 
the decay rates in the domain of outer communication, including its boundary T-L^, are 
in agreement with those obtained in [51 El H El El [12] • In particular, the asymptotic 
power law decay ~ T^" has been confirmed with n = /' + / + 3 for /> /' and with 
n = /' + / + 1 for / < /' as proposed in [5]. 

Due to the use of conformal compactification, the decay exponents could also be 
determined at J^"*" where our results support the analytically obtained predictions of 
[1| for the power low decay with n = / + 2 for / > /' and with n = /' for / < /' — 2. In 
particular, we have found that for the excitations /' = 4 and /' = 5 two of the lowest 
allowed excited modes, / = 0, 2 and / = 1, 3, have coinciding decay rates, n = 4 and 
n = 5, at , respectively. 

We extended the scope of numerical investigations by admitting not only axially 
symmetric configurations, with azimuthal index m = 0, but for two types of initial 
data specifications the cases with m = ±1 were also studied. Our findings, on the one 
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hand, confirmed the results of [5] indicating that the decay exponents n = n{l', I) are 
really independent of the azimuthal index. On the other hand, they provided some 
new insight about the nature of angular momentum transfers. 

There were certain situations where the method we applied was not providing 
definite values for the decay exponents. This happened either because the numerical 
accuracy was not sufficiently high or in some few cases because the local power index 
showed an unexpected location dependence and 'splitting' at the middle of the domain 
of outer communication. In particular, we found that in the latter cases the decay 
rates were smaller in the neighborhood of the black hole than in the distant region. 
We also found that this type of splitting sometimes remains and in some other cases 
disappears if non-axisymmetric (i.e. m = ±1) fields are taken into account. These 
cases definitely deserve further investigations. 

In addition to the fall off properties of the scalar field itself, we also investigated 
the energy and angular momentum conservation and transfers. 

Regarding the conservations, our numerical setup was found to be sufficiently 
accurate and even the convergence rate of our numerical code was found, in the 
initial burst phase, to be larger than four in accordance with the order of the applied 
finite differencing scheme. 

The transports were somewhat more interesting. We found that — in the tail 
phase — the decay exponents for the energy and angular momentum losses at J^~^ 
are smaller than at the horizon which means that considerably larger amount of en- 
ergy and angular momentum is radiated toward future null infinity than falling into 
the black hole. 

We also found that there is significant difference between the behavior of the 
energy and angular momentum transfers in the initial phase. As it is indicated by 
Figure [71 considerably larger part of the angular momentum of the scalar field gets 
into the black hole than out to J^^. What is even more interesting is that while at the 
beginning the inward moving impulse possesses almost the entire of the initial angular 
momentum, it starts to decrease as the impulse falls toward the black hole, and in the 
meanwhile, the role of the outward moving impulse gets to be more significant. As a 
result, at the beginning of the quasi-normal ringing phase considerably more angular 
momentum than energy stays in the domain of outer communication. 

It would be interesting to investigate the new features of the energy and angu- 
lar momentum transports in more detail. The above-mentioned ambiguities related 
to some of the fall off exponents would also deserve further investigations, which, 
however, are beyond the scope of the present paper as they require the use of much 
higher numerical accuracy and considerably longer time evolutions. It would also be 
important to derive analytic predictions for the late time decay exponents that we 
have found for stationary and for slowly falling off initial data. 
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Appendix A 



The coefficients appearing in the field equation (j8]) are as follows: 



arpR = 2{&MR'^ + 2a'^R^ + Qa^MR^ + AMR + SM^R^ + 2M - AMR^ - + SM'^R 

-8MR^ - AR^ - IQM^R^ - Ga'^MR'^ + 2a^MR^ - 2a^M - a'^R'^){l + R)R 
ar^ = AaR{-l + 2MR^ -2M){1 + R){1 + R^) 
aR^ = a{l + R^){1 + R){1 - 2R^ + R^) 

ar = — ^ (a^ + 2a'^M -2M + ia^ + 2a^M - 2M)R + (-4 + 5a^ + l2a^M - IAM)R^ 
1 + il^ 

+(-48M2 + l2a^M + ha^ -A- 30M)R'^ + (-Sa^ + 4 - 2Aa^M - ASM"^ + 10M)R^ 
+{A - Sa^ - 2Aa^M + 32M'^ + 26M)R^ + (-a^ + 32M^ + Aa'^M + 6M)R^ 
+(6M + Aa'^M + 16M2 - a^)/?^ + (IGM^ + Qa^M)R^ + Ga^MR^) 

2/1(1 + R ) 

+(12M - 12 + lOa^)/?"^ + (-12 + lOa^ - 24M)i?^ + (8 - 24M - 8a^)R^ 
+{AM - Sa^ + 8)R'^ + {AM - + 4)i?^ + (-a^ + 4 + 6M)ii9 
+(20^ + 6M)i?^° + 2a2i?") 

a(l+i?-i?2_^3)(i + ^2)2 

= R 

^ (a^ + (a^ - 2M)i? - 2MR^ - 2MR^ - {2a^ + 2M)R^ + {-2a^ + 2M)R^ 



ao 



2R^ 

+2MR^ + 2MR' + (2Af + a^)R^ + a^R^) 
OA = {I + R){1 + 2R^ + R^) 
aPj, = AMR + (32M2 - IQa^M + 20M - IGM^a^ + 4 - 4a2)i22 

+(64M^ - 16a^M + 44M + 4 + OGM^ - IGM^a^ - 40^)72^ 
+(64M^ + 32M2 + IGa^M + 2,2M'^a'^ - AM)R^ 

+(32M2a2 - 32M2 + IGa^M - QAM^)R^ + [-QAM^ - 16M'^a'^)R^ - IQM'^a^R' 

+-{a^ + a'^R + 2a^R^ + 20^/?^ + ^2^4 _^ ^2^5^) 
3 

a^^l = ^ -[a^ + a'^R + 2a^R^ + 2a^R^ + a^R^ + a^R^) 
3 V 5 
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Appendix B 



In this appendix, exphcit formulas are given for the integrands in the energy and 
angular momentum integrals J^rinE'^ and j^n^M^^, where f/ is a T = const or an 
R = const hyper surf ace. 

On a T = const hypersurface the integral j^j n^E^ can be written more explicitly 

as ^ 

di?d^d^ ^J\detgf\ (B.l) 

where g^^ denotes the metric on the 3-dimensional T = const surface obtained by 
the restriction of the full metric. The integral Jjj n^M^ on a T = const hypersurface 
can be written as 

di?d^d^ ^\detgf\ (B.2) 

Similarly, on an i? = const hypersurface, the the energy and the angular momentum 
integrals can be written as 

j dTd^d^ \Ad^ j dTdOdv \l\deigf\ (B.3) 

Let us introduce the quantities 

Gi - - (B.4) 

and 

ctr = 2{2a'^MR' + SM^/?^ + (6M - - 6a^M)R^ + {-AM - IGM^)/?^ 
+ (-8M - 4 + + Qa^M)R^ + {AM + m'^)R^ 

+ {-a^ + 2M - 2a^M)R)/{l + R^){1 + Rf (B.5) 
= A{-1 + 2MR^ -2M)aR/{l + Rf, cr^ = a (B.6) 



c^^ = {-2a^R^ -AMR + AM R^ + + AR^ + a^R^)/ 2{1 + R^) (B.7) 



di = -2{-16R'^M^a^ + {-WM^a^ -6AM^)R^ 
+ {-6AM^ + 32M^a^ - 32M^ + lQa^M)R^ 
+ (32M2a^ + UM^ + IQa^M + 32M2 - 4M)i?^ 
+ {-lQM^a^ - Aa^ + 44M - IGa^M + 96M^ + 4 + 6AM^)R^ 
+ {A - IGM^a^ + 20M + 32M2 - Aa^ - 16a^M)R'^ 
+AMR)/{1 + R)^{1 + R^) (B.8) 
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d2 = -2(a2i?5 + a^R^ + 2a^R' + 2a^R^ + a^R + a^)/{l + + i?^) (B.9) 

CTT = id, + ^)2V^Y,' + ^^d,Yl (B.IO) 
where yj* and 1^2° are the spherical harmonics 

The four integrands above can be expressed as 

V|det#|-^ = '^{Glcrrdr^dr^ 

-crr{G2^ + Gi{l - R)dR^){G2^ + Gi{l - R)dR^) 
-CR^Giil - i?)a^^(G2* + Gi{l - R)dR^) 
-CR^Giil - R)d^^{G2^ + - R)dR^) 
+(1 + i?2)G?(^As2^ + ^As2^)/(1 + (B.12) 



+CTRG,d^'^{-^^ + GiOr^) 
+CTRGid^^{-^'^ + GiOr^) 

i — K 

+2cT^Gld^^d^'^) (B.13) 



V'ld^t^ = ^ (cRRG^il - R)dT^{G2^ + - R)dR^) 

+c^^G'i(l - i?)9T^(G2^ + Gi(l - i?)^^^) 

+cr^G\{\ - Rfid^dr"^ + (B.14) 
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sin 6 



4 



+crrGi{1 - R)d^^{G2^ + G,{l - R)dn^) 
+2cR^Gl{l - Rfd^^d^'^) 



(B.15) 



The integrals over 6 G [0, vr] and ip G [0, 2tt] can be evaluated by considering that 
Jq d9 Jq^ dip sinO f{6,(p) = /g, where f{0,ip) is an arbitrary function and /q is the 
coefficient of the / = m = mode in its spherical harmonic expansion. The integrals 
over T and R are evaluated numerically. 
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